%% create discretized maps and graphs of road network
tic
clear
close all
clc
set(0,'DefaultFigureWindowStyle','docked','DefaultFigureVisible','on')

% -----------------------
% DEFINE FOLDER LOCATIONS
folders;

%% code control

    % which cities to keep?
    popthreshold = 0.1;       % only load places with population above threshold
    pick_largest = 0.6;       % =1 to use most populated place in a cell to locate the node in the cell, otherwise choose centroid
    max_N_cities = inf;       % keep the max_N_cities largest cities
    N_cat = 10;                % for plotting only. number of population quantiles is N_pop_cat+1.

    % discretization
    % these parameters control the discretization of the road network
    diagonals = 1;    % include diagonals 

    % keep largest connected component?
    largest_connected = 1;
    
    % testing income data 
    winsorize = 1;

    % when to build for forests/reserve areas 
    build_thresh = 80; 
    default_cellsize = 0.5;   % default size is 50km by 50km (0.5)
    adjust_cellsize = 1; % adjust cell size based on country size 
    
    % add unconnected nodes?
    connect_all_nodes = 1;   % number of times the original network is made denser (times unconnected nodes are connected)
    tol = 0.0001; 
    save_actual = 0;

    % which pop source to use?
    WP = 0; % 0 - default, 1 - WorldPop
    
    % which measure of infrastructure quality?
    osm = 1; % 0 - default, 1 - OSM
    
    % stock all code control variables
    code_control.build_thresh = build_thresh;
    code_control.datafolder_ne = datafolder_ne;
    code_control.popthreshold = popthreshold;       
    code_control.pick_largest = pick_largest;        
    code_control.max_N_cities = max_N_cities;    
    code_control.N_cat = N_cat;               
    code_control.diagonals = diagonals;
    code_control.winsorize = winsorize;

%% load country list
country_list_short;


%% loop through selected countries
for country_n= 1:Ncountries
 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% load international boundaries %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
    country_icc = char( countries(country_n) );  % country ICC code
    country = icc2name( country_icc );           % country name
    
    % get threshold from 
    x_ver_hor = get_threshold( country );
    x_diag  = x_ver_hor;
    
    % add parameters to code control 
    code_control.x_ver_hor = x_ver_hor;  
    code_control.x_diag = x_diag;
    
    if ~ispc()   %% we are in a mac

        % folder to save the outcome
        datafolder_LAC_country =  [ datafolder_LAC,country_icc,'/'];

    else   %% we are in a pc   

        % folder to save the outcome
        datafolder_LAC_country =  [ datafolder_LAC,country_icc,'\' ];

    end
    
    % load country boundaries
    disp( ['Country: ',country] )
    file  = [ datafolder_LAC_country,'country_bounds.mat'];
    load(file)
           
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% load populated places %%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   % load Population Data
   file  = [ datafolder_LAC_country,'popplaces_sedac.mat'];
   temp = load(file);   
   places_sedac = temp.places_sedac;
   clear temp
     
   % total population
   totpop = sum( cell2mat( {places_sedac.population} ) );
   
    % World Pop
    file  = [ datafolder_LAC_country,'popplaces_worldpop.mat'];
    load( file,'places_WorldPop' );
    
    % LandScan
    file  = [ datafolder_LAC_country,'popplaces_landscan.mat'];
    load( file,'places_LandScan' );

%%  %%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% bring income data %%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    
    file  = [ datafolder_LAC_country,'places_gecon.mat'];
    load( file,'places_gecon' );
    
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% bring altitude data %%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    file  = [ datafolder_LAC_country,'relief_etopo.mat'];
    load( file,'relief_etopo' );      
       
%%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%% discretize the map %%%% 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    % places_grid is a structure with points
        % with the same structure as places, but constructeed on averages
        % over the discrete grid. It also includes neighbors of each point.
        % gridmap is a structure with polygons
   
    if strcmp(country, 'Peru')
        
        % load in forest data 
        file  = [ datafolder_LAC_country,'gfc.mat'];
        load( file,'gfc' );     
        [ code_control,places_grid,gridmap,edges,places_gecon ] = ...
            create_gridmap_forest( code_control,country_bounds, gfc, places_sedac,places_WorldPop, places_LandScan,relief_etopo,places_gecon,default_cellsize,adjust_cellsize );

    elseif strcmp(country, 'Brazil')

            % brazil's gridmap is different 
            [ code_control,places_grid,gridmap,edges,places_gecon ] = ...
            brazil_grid( code_control,country_bounds,places_sedac,relief_etopo,places_gecon);
    else

        % standard grid creation 
        [ code_control,places_grid,gridmap,edges,places_gecon ] = ...
            create_gridmap( code_control,country_bounds,places_sedac, ...
            places_WorldPop, places_LandScan, ...
            relief_etopo,places_gecon,default_cellsize,adjust_cellsize, country );
        
        % we allow building everywhere 
        for j=1:length(gridmap)
            gridmap(j).build = 1;
        end
        
    end
    
    % Need to remove area in natural reserve areas for Chile
    if strcmp(country, 'Chile')
        
        % load in the nature preserve shapefile
        marine = @(code) test_marine(code);
        
        % load natural reserve areas 
        chile_nature_preserve = shaperead([path_raw_mapdata, '/nature_shp/WDPA_Apr2020_CHL/WDPA_Apr2020_CHL_shapefile_polygons'], 'Selector',{marine,'MARINE'}); 

        % first simplify the shapefile 
        % for each cell, compute share of cell that is a nature preserve
        nature_shp_X = cell2mat({chile_nature_preserve.X});
        nature_shp_Y = cell2mat({chile_nature_preserve.Y});

        % reduce detail in shapefile 
        [nature_Y1, nature_X1, cerr, tol] = reducem( nature_shp_Y', nature_shp_X', 1.6548e-04);

        % check if the center of the cell is in the shapefile 
        polyarray2 = polyshape(nature_X1, nature_Y1);

        in_shp = zeros( length(gridmap), 1);
        for j=1:length(gridmap)

            [ j length(gridmap) ]

            polyarray1 = polyshape([gridmap(j).X], [gridmap(j).Y]);
            poly1 = [polyarray1 polyarray2];
            polyout = intersect(poly1);
            in_shp(j) = area(polyout)/area(polyarray1);
            gridmap(j).build = in_shp(j)<0.80;
            places_grid(j).build = gridmap(j).build;
        end
    
    end
       
   % set population and income per capita depending on pop source selected
   if WP == 1 
        for j=1:length(places_grid)
            places_grid(j).population = gridmap(j).population_WP;
            places_grid(j).y = gridmap(j).yWP;
        end        
   end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%% load in travel speeds %%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if osm == 1 
        file  = [ datafolder_OSM, country_icc, 'coordinates_wtimes.csv'];
        osm_speeds = readtable(file);
        for n=1:height(osm_speeds)
            if strcmp(osm_speeds.speed(n), 'NA') || isnan(osm_speeds.speed(n))
                osm_speeds.speed(n) = 0;
            end
        end
    
        speeds  = cell2mat({osm_speeds.speed})';
        
        % Missing links that are neighbors can always be walked between
        % at a rate of 4km/hr
        speeds(speeds<4) = 4;
        speeds = speeds';
    end
        
%%  %%%%%%%%%%%%%%%%%%%%
    %%%% load roads %%%% 
    %%%%%%%%%%%%%%%%%%%%
    file = [datafolder_LAC_country, country_icc, '_GRIP.shp'];
    roads_temp1 = shaperead( file);     

    % export roads to mat file for plotting, summary statistics 
    if save_actual == 1
        roads_export = assign_weights_to_roads_LAC( roads_temp1 );
        save( [ path_save_grids,country,'_roads.mat'], 'roads_export', '-v7.3');  
    end

    % bound roads - this is important, as we learned...
    use_hull = 0;
    roads_temp2 = bound_roads( roads_temp1,country_bounds,use_hull );
    
    % add weight to each segment of road based on its characteristics to compute shortest paths
    roads = assign_weights_to_roads_LAC( roads_temp2 );

    clear roads_temp1;
    clear roads_temp2;
    clear roads_export;

    %% transform the roads into graph
    [ road_network,unique_nodes,unique_nodes_nat ] = transform_network_to_graph( code_control,roads,country );

        %% Make the network denser to solve the connectivity problem
    if connect_all_nodes==1 

        if strcmp(country, 'Brazil')
            tol = tol*10; % Brazil's road network has lots of gaps 
        end

        road_network = denser_network(road_network, unique_nodes, tol, min(road_network.Edges.Weight));
    end

    %% Keep the largest connected component of road network
    if largest_connected==1   
        
        fprintf('Keeping largest connected component')
        connected_component = conncomp( road_network);
        idx = connected_component == mode(connected_component);
        connected_roads = subgraph(road_network, idx);
        unique_nodes = unique_nodes( idx,: ); 

        % get the national component   
        % get the indicies of unique_nodes_nat that appear in the 
        unique_nodes_nat = intersect(cell2mat({unique_nodes_nat}), cell2mat({unique_nodes}), 'rows');
        road_network2 = connected_roads;
    
    else 
        
        road_network2 = road_network;
    end        
        
    % move cities and discretize roads     
    places2 = move_cities( gridmap,places_grid,unique_nodes,unique_nodes_nat);  
    
    % no diagonals here since it is not a true grid 
    if strcmp(country, 'Brazil')

       code_control.diagonals = 0; 
    
    end
    
    % run code to discretize the road network 
    if osm == 0

        [ discretized_roads,graph_export,unique_edges, dist_info ] = discretize_roads( code_control,gridmap,places2,road_network2,edges,unique_nodes, country );        
   
    elseif osm== 1

        [ discretized_roads,graph_export,unique_edges, dist_info ] = discretize_roads_speed( code_control,gridmap,places2,road_network2,edges,unique_nodes, country, speeds );        

    end
    
    country_graph.places_grid=places_grid;
    country_graph.country_bounds=country_bounds;
    country_graph.gridmap=gridmap;
    country_graph.places2=places2;
    country_graph.discretized_roads = discretized_roads;
    country_graph.graph_export = graph_export;
    country_graph.edges = edges;  
    country_graph.unique_edges = unique_edges;      
    country_graph.code_control = code_control;
    country_graph.unique_nodes = unique_nodes; 
    country_graph.dist_info = dist_info;

    % export graph for use in calibration 
      save( [ path_save_grids,country,'_grid_',...
        num2str( x_diag ),'_',...
        num2str( x_ver_hor ), '_',...
        num2str( default_cellsize ),'_connected', ...
        num2str( largest_connected ), ...
        '_WP', num2str(WP), '_osm', num2str(osm),'.mat' ],'country_graph' );   

end



        